1 Introduction

<<<<<<< HEAD

=======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

1.1 Why do we want to assess community changes ?

Global changes have affect a tremendous effect communities
Understanding assembly and disassembly rules
Give a more complex and complete picture of what is going on ?

<<<<<<< HEAD

=======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

1.2 Documenting community changes

Many papers report temporal trends of local species richness and turnover
- Abundance (Klink et al. 2020)
- Species richness and or turnover: Blowes et al. (2019)
- Species richness and or turnover: (Vellend et al. 2013, @dornelas_assemblage_2014, @blowes_geography_2019)
- They describe basically no trends in species richness but high turnover
- Species turnover vary a lot according to geography and taxons (Blowes et al. 2019)
- Species turnover (as opposite of nestedness (baselga_temporal_2013?)) has =======
- Species turnover (as opposite of nestedness (???)) has >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432 been reported to be the major change, i.e. change in species identity, which is coherent with the absence of species richness trends

<<<<<<< HEAD

=======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

1.3 Research gap

Biodiversity trends are all about geographical variation, i.e. few ecological drivers are involved

<<<<<<< HEAD

1.4 Research questions

=======

1.4 Research questions

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

1.4.1 Drivers of temporal changes

What are the drivers of community changes ?
- Climate change: more appearance and disappearance in the upper basin
- More species turnover far from the source because of dispersal ability

cap_turn <- paste0(
  "The hypothesis is that sites located at the middle of the stream network
  will experience the more composition change through time because of two
  opposites constraints. Climate change is expected to drive species upward in
  the stream network but an opposite constrain gathered into hydrological
  contrain, such as water flowand obstacles to migration will restrain the
  upward migration. The sum of those constrains is hypothesized to result in a
  quadratic relationship between jaccard turnover and distance from source.")

<<<<<<< HEAD

knitr::include_graphics(here("doc/fig/turnover_hyp.pdf"))

The hypothesis is that sites located at the middle of the stream network
=======
<pre class=

knitr::include_graphics("fig/turnover_hyp.png")

The hypothesis is that sites located at the middle of the stream network
>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
will experience the more composition change through time because of two
opposites constraints. Climate change is expected to drive species upward in
the stream network but an opposite constrain gathered into hydrological
contrain, such as water flowand obstacles to migration will restrain the
upward migration. The sum of those constrains is hypothesized to result in a
quadratic relationship between jaccard turnover and distance from source.

Figure 1.1: The hypothesis is that sites located at the middle of the stream network will experience the more composition change through time because of two opposites constraints. Climate change is expected to drive species upward in the stream network but an opposite constrain gathered into hydrological contrain, such as water flowand obstacles to migration will restrain the upward migration. The sum of those constrains is hypothesized to result in a quadratic relationship between jaccard turnover and distance from source.

<<<<<<< HEAD

=======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

1.4.2 Within and basin scale

Temporal homogeneization at the basin scale?

dis_cap <- paste0("If sites within basin becomes more dissimilar from their
  reference (first year of observation), we
  should observe at the same time an homogeneization of sites at the basin scale")

ggplot(tibble(x = c(0, .1)), aes(x)) +
  stat_function(fun = function(x) x, geom = "line") +
  labs(
    x = "Average temporal dissimalitiry within basin",
    y = "Temporal similarity at the basin scale"
  )

<<<<<<< HEAD

If sites within basin becomes more dissimilar from their
=======
<img src=

Figure 1.2: If sites within basin becomes more dissimilar from their reference (first year of observation), we should observe at the same time an homogeneization of sites at the basin scale

Product:
- Map of temperature and distance to effect on turnover by basin ?

<<<<<<< HEAD

=======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

1.5 Fish freshwater in streams

Streams are a good model to study turnover:
- They are spatially delimited at the basin level
- Stream fishes are strongly constrained by hydrological constrains, which are highly variable at a small spatial scale, making possible to ecological drivers at multiple places, here the basin scale.

<<<<<<< HEAD

=======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

2 Methods

2.1 Data selection

Unique protocol and unique unit of abundance by site
One sampling by year and site
In each site, samplings should be in the median (more or less one month and half, i.e. 1.5 month)
5 samplings by site, span 10 years

loc <- filtered_dataset$location %>%
  left_join(filtered_dataset$site_quali, by = "siteid") %>%
  st_as_sf(coords = c("longitude", "latitude"),
  crs = 4326)

fig_sites <- paste0(
  "Australia absence because span is superior to 10 years of 100s of sites but number of samplings < 10.")

world <- ne_countries(scale = "medium", returnclass = "sf") %>% 
  st_transform(crs = 4326)
bb <- st_bbox(loc)
ggplot(data = world) +
    geom_sf() +
    geom_sf(data = loc, aes(color = protocol), shape = 1) +
#    coord_sf(xlim = bb[c("xmin", "xmax")],
#      ylim = bb[c("ymin", "ymax")], 
#      datum = 4326, expand = TRUE) +
    theme(legend.position = "bottom") +
    labs(title = paste0("Number of sites: ", nrow(loc)))

<<<<<<< HEAD

Australia absence because span is superior to 10 years of 100s of sites but number of samplings < 10.

=======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

Figure 2.1: Australia absence because span is superior to 10 years of 100s of sites but number of samplings < 10.

2.2 Biodiversity

2.2.1 Quantity

Hill numbers
- $H^0$: Species richness
- $H^1$: Shannon
- $H^2$: Simpson

We used also hill numbers with coverage-based correction proposed by (chao_coverage-based_2014?)

Evenness:
- Pielou: $J = H^1 / log(S)$

skim(analysis_dataset, starts_with("chao_"), "species_nb", "shannon", "simpson",
"evenness")

Table 2.1: Data summary

Name	analysis_dataset
Number of rows	62174
Number of columns	57
_______________________
Column type frequency:
numeric	8
________________________
Group variables	None

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
chao_richness	1	5.12	5.45	1	1.94	3.14	6.44	85.71	<<<<<<< HEAD <U+2587><U+2581><U+2581><U+2581><U+2581> ======= ▇▁▁▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
chao_shannon	1	2.79	2.59	0	1.35	2.13	3.63	33.63	<<<<<<< HEAD <U+2587><U+2581><U+2581><U+2581><U+2581> ======= ▇▁▁▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
chao_simpson	1	2.23	1.92	0	1.20	1.79	2.89	40.00	<<<<<<< HEAD <U+2587><U+2581><U+2581><U+2581><U+2581> ======= ▇▁▁▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
chao_evenness	1	1.93	1.54	0	1.47	2.06	2.47	38.25	<<<<<<< HEAD <U+2587><U+2581><U+2581><U+2581><U+2581> ======= ▇▁▁▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
species_nb	1	5.60	5.92	1	2.00	3.00	7.00	72.00	<<<<<<< HEAD <U+2587><U+2581><U+2581><U+2581><U+2581> ======= ▇▁▁▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
shannon	1	0.86	0.65	0	0.34	0.76	1.29	3.54	<<<<<<< HEAD <U+2587><U+2586><U+2583><U+2581><U+2581> ======= ▇▆▃▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
simpson	1	0.42	0.28	0	0.17	0.48	0.66	0.96	<<<<<<< HEAD <U+2587><U+2585><U+2587><U+2587><U+2583> ======= ▇▅▇▇▃ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
evenness	1	0.54	0.30	0	0.33	0.63	0.77	1.00	<<<<<<< HEAD <U+2585><U+2582><U+2585><U+2587><U+2585> ======= ▅▂▅▇▅ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

2.2.2 Betadiversity

Betadiversity metrics compares the composition of two communities. In our case, we compare the composition of each time step relatively to the baseline, i.e. the first year of sampling.

Species exchange index: (Hillebrand et al. 2018)
- Presence-absence: $J = SER_r = \dfrac{S_{imm} + S_{lost}}{S_{tot}}$
- Relative abundance: $SER_a = \dfrac{\sum_i (p_i - p^{\prime}_i)^2}{\sum_i p_i^2 + \sum_i p^{\prime2}_i - \sum_i p_i p^{\prime}_i}$
- $S_{imm}$, $S_{lost}$, $S_{tot}$: number of immigrant, lost and total species respectively
- $i$: species $i$, $p$: relative abundance, $\prime$: second community
- $J$: Jaccard index

Turnover metrics based on presence/absence can be decomposed in appearance/disappearance and nestedness / turnover.

appearance: $\dfrac{S_{imm}}{S_{tot}}$, disappearance: $\dfrac{S_{lost}}{S_{tot}}$
Turnover: $J_t = \dfrac{2 * min(S_{lost}, S_{imm})}{S_{common} + (2 * min(S_{lost}, S_{imm})}$
Nestedness: $SER_r - J_t$

skim(analysis_dataset, "jaccard", "turnover", "nestedness", "appearance",
  "disappearance", "hillebrand")

Table 2.2: Data summary

Name	analysis_dataset
Number of rows	62174
Number of columns	57
_______________________
Column type frequency:
numeric	6
________________________
Group variables	None

Variable type: numeric

skim_variable	complete_rate	mean	sd	p25	p50	p75	p100	hist
jaccard	1	0.65	0.27	0.50	0.67	1.00	1.00	<<<<<<< HEAD <U+2581><U+2583><U+2587><U+2586><U+2587> ======= ▁▃▇▆▇ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
turnover	1	0.16	0.25	0.00	0.00	0.33	1.00	<<<<<<< HEAD <U+2587><U+2582><U+2581><U+2581><U+2581> ======= ▇▂▁▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
nestedness	1	0.19	0.22	0.00	0.10	0.33	0.96	<<<<<<< HEAD <U+2587><U+2582><U+2582><U+2581><U+2581> ======= ▇▂▂▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
appearance	1	0.14	0.16	0.00	0.10	0.24	0.92	<<<<<<< HEAD <U+2587><U+2583><U+2581><U+2581><U+2581> ======= ▇▃▁▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
disappearance	1	0.11	0.14	0.00	0.00	0.17	0.86	<<<<<<< HEAD <U+2587><U+2582><U+2581><U+2581><U+2581> ======= ▇▂▁▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
hillebrand	1	0.68	0.33	0.42	0.79	0.99	1.00	<<<<<<< HEAD <U+2582><U+2582><U+2582><U+2582><U+2587> ======= ▂▂▂▂▇ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

<<<<<<< HEAD

3 Environment

=======

3 Environment

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

3.1

tar_load(riveratlas_site)
riv <- riveratlas_site %>%
  select(all_of(c("siteid", setNames(get_river_atlas_significant_var(), NULL)))) %>%
  mutate(across(starts_with("tmp_"), ~.x / 10)) %>%
  st_drop_geometry()

riv %>%
  rename(get_river_atlas_significant_var()) %>%
  pivot_longer(cols = -siteid, names_to = "variable", values_to = "values") %>%
  ggplot(aes(x = values)) +
  geom_histogram() +
  facet_wrap(~variable, scales = "free")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

<<<<<<< HEAD

=======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

4 Statistical analysis

skim(modelling_data)

Table 4.1: Data summary

Name	modelling_data
Number of rows	54626
Number of columns	<<<<<<< HEAD 20 ======= 18 >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
_______________________
Column type frequency:
character	1
factor	1
numeric	<<<<<<< HEAD 18 ======= 16 >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
siteid	0	1	2	6	0	4916	0

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
main_bas	0	1	FALSE	291	708: 2696, 208: 2400, 208: 2252, 208: 1566

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
year	1	2004.86	8.57	1955.0	1999.00	2006.00	2012.00	2019.00	<<<<<<< HEAD <U+2581><U+2581><U+2582><U+2587><U+2587> ======= ▁▁▂▇▇ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
year_nb	1	11.17	8.66	0.0	4.00	10.00	17.00	63.00	<<<<<<< HEAD <U+2587><U+2583><U+2581><U+2581><U+2581> ======= ▇▃▁▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
scaled_dist_up_km	1	0.00	1.00	-0.4	-0.35	-0.28	-0.11	13.99	<<<<<<< HEAD <U+2587><U+2581><U+2581><U+2581><U+2581> ======= ▇▁▁▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
span	1	22.45	8.46	10.0	16.00	22.00	28.00	64.00	<<<<<<< HEAD <U+2587><U+2587><U+2582><U+2581><U+2581> ======= ▇▇▂▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
jaccard_scaled	1	0.65	0.27	0.0	0.50	0.66	1.00	1.00	<<<<<<< HEAD <U+2581><U+2585><U+2587><U+2586><U+2587>
jaccard_dis	1	0.35	0.27	0.0	0.00	0.34	0.50	1.00	<U+2587><U+2586><U+2587><U+2583><U+2581> ======= ▁▅▇▆▇ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
turnover	1	0.16	0.25	0.0	0.00	0.00	0.33	1.00	<<<<<<< HEAD <U+2587><U+2582><U+2581><U+2581><U+2581> ======= ▇▂▁▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
nestedness	1	0.19	0.22	0.0	0.00	0.10	0.33	0.96	<<<<<<< HEAD <U+2587><U+2582><U+2582><U+2581><U+2581> ======= ▇▂▂▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
species_nb	1	5.88	6.05	1.0	2.00	4.00	7.00	72.00	<<<<<<< HEAD <U+2587><U+2581><U+2581><U+2581><U+2581> ======= ▇▁▁▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
log_species_nb	1	1.38	0.87	0.0	0.69	1.39	1.95	4.28	<<<<<<< HEAD <U+2587><U+2587><U+2586><U+2582><U+2581> ======= ▇▇▆▂▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
chao_richness	1	5.36	5.56	1.0	1.96	3.37	6.90	85.71	<<<<<<< HEAD <U+2587><U+2581><U+2581><U+2581><U+2581> ======= ▇▁▁▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
hillebrand	1	0.67	0.33	0.0	0.41	0.78	0.99	1.00	<<<<<<< HEAD <U+2582><U+2582><U+2582><U+2582><U+2587>
jaccard_dis_scaled	1	0.35	0.27	0.0	0.00	0.34	0.50	1.00	<U+2587><U+2586><U+2587><U+2583><U+2581> ======= ▂▂▂▂▇ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
nestedness_scaled	1	0.19	0.22	0.0	0.00	0.10	0.33	0.96	<<<<<<< HEAD <U+2587><U+2582><U+2582><U+2581><U+2581> ======= ▇▂▂▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
turnover_scaled	1	0.16	0.25	0.0	0.00	0.00	0.33	1.00	<<<<<<< HEAD <U+2587><U+2582><U+2581><U+2581><U+2581> ======= ▇▂▁▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
hillebrand_scaled	1	0.67	0.33	0.0	0.41	0.78	0.99	1.00	<<<<<<< HEAD <U+2582><U+2582><U+2582><U+2582><U+2587> ======= ▂▂▂▂▇ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
log1_year_nb	1	2.15	0.96	0.0	1.61	2.40	2.89	4.16	<<<<<<< HEAD <U+2583><U+2583><U+2587><U+2587><U+2581> ======= ▃▃▇▇▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
one	1	1.00	0.00	1.0	1.00	1.00	1.00	1.00	<<<<<<< HEAD <U+2581><U+2581><U+2587><U+2581><U+2581> ======= ▁▁▇▁▁ >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

<<<<<<< HEAD

=======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

4.1 Species composition change

In each site and at each sampling point, we computed the dissimilarity (Jaccard, Appearance, Disappearance, Turnover & Nestedness) between the time step and the reference time, i.e. the first year of sampling. Regarding the reference time step, there are two different strategies. First is to include it, i.e. setting a dissimilarity of 0 at each beginning of a time series (Blowes et al. 2019). Another strategy is to exclude the first year from the analysis (Dornelas et al. 2014). The first option logically constrains the analysis toward increasing dissimilarity as the second can capture recovery from perturbation soon after the first samplings. I guess that the first option makes more sense more it removes to sensibility to community changes happening just after the references year.

<<<<<<< HEAD =======

4.2 Model

Simple lm by site
Model all in one:

4.2.1 Turnover

Jaccard similarity, nestedness and turnover were modelled in two ways, Beta and gaussian regression. Beta regressions are more appropriate for distribution borned by [0,1] and are linearized by a logit link. Logit links have the cons that coefficients are less easy to interpret. So we run both gaussian and beta regressions and investigate the relationship among coefficients of those type modelling types. Values of turnover were slightly transformed such as values of one and zero to fit beta regression requirements ($x' = \dfrac{x \times (N - 1) + 0.5}{N}$).

The model is defined like this:

paste0(
  "jaccard_scaled ~
  0 + year_nb * scaled_dist_up_km +
  (0 + year_nb | main_bas/siteid) +
  (0 + year_nb | span) +
  (0 + year_nb:scaled_dist_up_km | main_bas)"
)
#> [1] "jaccard_scaled ~\n  0 + year_nb * scaled_dist_up_km +\n  (0 + year_nb | main_bas/siteid) +\n  (0 + year_nb | span) +\n  (0 + year_nb:scaled_dist_up_km | main_bas)"

year_nb: number of year since the beginning of the monitoring (site scale)
scaled_dist_up_km: Distance from the source (centered and scaled)
main_bas: Hydrographic basin
siteid: site
span: Number of years between the end and the beginning of the monitoring

Our objective was to assess the temporal turnover and its ecological drivers. We set intercept to equal 1.0 because year_nb = 0 at the beginning of each series, jaccard is always equal to 1.0. We estimated the effects of the basin and of the site on the temporal trends (0 + year_nb | main_bas/siteid), the site effect being nested in the basin effect. We also evaluated the effects of span on the temporal trends, because it is generally expected that longer timeseries displays stronger temporal trends. Finally, we acknowledged that the effect of the distance from source can be different among basin.

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

4.2 Model

Simple lm by site
Model all in one:

4.2.1 Turnover

The model is defined like this:

paste0(
  "jaccard_scaled ~
  0 + year_nb * scaled_dist_up_km +
  (0 + year_nb | main_bas/siteid) +
  (0 + year_nb | span) +
  (0 + year_nb:scaled_dist_up_km | main_bas)"
)
#> [1] "jaccard_scaled ~\n  0 + year_nb * scaled_dist_up_km +\n  (0 + year_nb | main_bas/siteid) +\n  (0 + year_nb | span) +\n  (0 + year_nb:scaled_dist_up_km | main_bas)"

year_nb: number of year since the beginning of the monitoring (site scale)
scaled_dist_up_km: Distance from the source (centered and scaled)
main_bas: Hydrographic basin
siteid: site
span: Number of years between the end and the beginning of the monitoring

5 Results

5.1 Direction of trends

tar_load(var_analysis)
ti <- map_dfr(rigal_trends, ~tabyl_df(x = .x, group = "direction"),
  .id = "response"
)
ti %>%
  filter(direction != "Total", response %in% var_analysis) %>%
  select(response, direction, percent) %>%
  mutate(response = str_replace_all(response, get_var_replacement())) %>%
  pivot_wider(names_from = "direction", values_from = "percent") %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "bordered", "hover"))

response	stable	increase	decrease
Species richness	81.1%	11.5%	5.7%
Log species richness	79.5%	11.1%	5.2%
Chao species richness	83.2%	9.4%	6.1%
SER_a (rel abundance)	73.5%	0.7%	24.5% <<<<<<< HEAD
Jaccard (binary, dissimilarity)	67.4%	26.8%	0.3% ======= >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
Nestedness (jaccard)	80.3%	13.1%	0.5%
Turnover (jaccard)	53.9%	12.2%	0.2%

ti <- map_dfr(rigal_trends, ~tabyl_df(x = .x, group = "shape_class"),
  .id = "response"
)

ti %>%
  filter(shape_class != "Total") %>%
  select(response, shape_class, percent) %>%
  mutate(response = str_replace_all(response, get_var_replacement())) %>%
  pivot_wider(names_from = "shape_class", values_from = "percent") %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "bordered", "hover"))

response	stable_constant	increase_constant	decrease_constant	stable_convex	stable_concave	decrease_decelerated	increase_accelerated	decrease_accelerated	increase_decelerated	NA_constant
Total abundance	72.7%	9.1%	6.2%	6.0%	4.1%	1.0%	0.8%	0.1%	0.1%	NA
Log Total Turnover (jaccard) (codyn) abundance	69.5%	10.4%	6.9%	4.8%	6.3%	0.3%	0.2%	0.4%	1.2%	0.0%
Species richness	70.7%	9.5%	4.8%	4.3%	6.0%	0.7%	0.8%	0.3%	1.1%	1.8%
Log species richness	69.3%	9.3%	4.4%	3.4%	6.9%	0.5%	0.7%	0.3%	1.1%	4.3%
Chao species richness	74.1%	7.8%	5.3%	4.6%	4.5%	0.7%	0.9%	0.1%	0.7%	1.3%
Chao Shannon	71.1%	8.5%	5.1%	4.2%	5.5%	0.4%	0.7%	0.2%	0.6%	3.7%
Chao Simpson	72.2%	7.8%	5.2%	4.3%	5.2%	0.5%	0.6%	0.1%	0.5%	3.7%
Chao Evenness	77.5%	5.5%	4.6%	3.3%	4.2%	0.3%	0.4%	0.1%	0.4%	3.7%
Jaccard (binary, similarity)	55.1%	0.3%	22.0%	14.3%	1.0%	4.5%	NA	0.7%	0.3%	1.8%
Horn (binary, similarity)	58.7%	0.3%	22.1%	11.9%	1.1%	3.1%	NA	0.8%	0.3%	1.8%
Chao (binary, similarity)	59.3%	0.9%	7.6%	9.2%	3.9%	3.7%	NA	1.3%	1.7%	12.4%
SER_a (rel abundance)	60.3%	0.5%	19.5%	11.8%	1.5%	3.9%	0.0%	1.2%	0.1%	1.3%
Total Turnover (jaccard) (codyn)	56.4%	22.1%	0.3%	0.8%	11.1%	NA	0.9%	NA	3.0%	5.5%
Appearance	58.9%	16.7%	0.4%	0.9%	8.5%	NA	1.2%	0.0%	1.8%	11.6%
disAppearance	53.0%	12.5%	0.2%	1.1%	7.1%	NA	0.6%	NA	1.7%	23.7%
Evenness	73.6%	6.7%	5.6%	4.1%	4.8%	0.4%	0.5%	0.2%	0.5%	3.6%
Shannon	70.3%	8.9%	5.6%	4.0%	5.7%	0.3%	0.6%	0.2%	0.7%	3.7%
Simpson	72.2%	7.8%	5.1%	3.7%	5.6%	0.3%	0.5%	0.3%	0.8%	3.6%
Jaccard (binary, dissimilarity)	53.0%	22.0%	0.3%	0.7%	13.7%	NA	0.7%	NA	4.2%	5.5%
Nestedness (jaccard)	71.7%	11.2%	0.5%	1.0%	7.6%	NA	0.8%	NA	1.1%	6.2%
Turnover (jaccard)	48.0%	10.5%	0.2%	1.0%	4.9%	NA	0.6%	0.0%	1.1%	33.7%

5.2 Average

slope <- map_dfr(rigal_trends,
  ~.x %>%
  select(siteid, linear_slope),
.id = "response"
)

rigal_trends_df <-  map2_dfr(
  rigal_trends, names(rigal_trends),
  ~.x %>% mutate(variable = .y)
  ) %>%
select(variable, siteid, linear_slope) %>%
pivot_wider(names_from = "variable", values_from = "linear_slope")

summary_slope <- slope %>%
  filter(response %in% var_analysis) %>%
  group_by(response) %>%
  summarise(summ = list(enframe(summary_distribution(linear_slope, na.rm = TRUE)))) %>%
  unnest(cols = summ) %>%
  pivot_wider(names_from = "name", values_from = "value") %>%
  select(response, mean, median, sd)
summary_slope %>%
  mutate(response = get_var_replacement()[response]) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "bordered", "hover")) %>%
  scroll_box(width = "100%", height = "1000px")

response	mean	median	sd
Chao species richness	0.0160926	0.0027879	0.1864729
SER_a (rel abundance)	-0.0184716	-0.0134346	0.0216062 <<<<<<< HEAD
Jaccard (binary, dissimilarity)	0.0180388	0.0152704	0.0171575 ======= >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
Log species richness	0.5690175	0.1305570	3.5608475
Nestedness (jaccard)	0.0090678	0.0053424	0.0154355
Species richness	0.0233072	0.0038939	0.1763604
Turnover (jaccard)	0.0089710	0.0022039	0.0161379

p_jt_sup_ne <- sum(abs(rigal_trends_df$turnover) > abs(rigal_trends_df$nestedness)) /
  nrow(rigal_trends_df)

To compare with Dornelas et al. (2014):
- Species richness:
  - Dornelas et al. (2014): 1557 measurements of species richness in two consecutive times, 629 (40%) increase, 624 (40%) decrease, and 304 (20%) do not change.
  - Our study: 13% increase, 7% decrease, 78% do not change, median: + 131% species by decade (i.e from 2 to 3 species, 10 to 15 species in 10 years)
- Jaccard:
  - 10% species replaced per decade (Dornelas et al. 2014)
  - 28% of species replaced per decade (Blowes et al. 2019)
  - our study: of species per decade
- Turnover / nestedness:
  - Blowes et al. (2019): 97% of study shows turnover exceeding nestedness
  - our study: 43%
- Jaccard vs dominance-based SER:
  - low correlation: $0.52$

tar_load(gaussian_jaccard_tmb)
library(glmmTMB)
confint(gaussian_jaccard_tmb$mod[[1]], "year_nb", component = "cond")
<<<<<<< HEAD
#>              2.5 %     97.5 %   Estimate
#> year_nb 0.02280811 0.02984302 0.02632557
summary(gaussian_jaccard_tmb$mod[[1]])
#>  Family: gaussian  ( identity )
#> Formula:          
#> jaccard_dis_scaled ~ 0 + year_nb/scaled_dist_up_km + (0 + year_nb |  
#>     main_bas/siteid) + (0 + year_nb | span) + (0 + year_nb:scaled_dist_up_km |  
#>     main_bas)
=======
#>             2.5 %      97.5 %    Estimate
#> year_nb -0.029843 -0.02280809 -0.02632554
summary(gaussian_jaccard_tmb$mod[[1]])
#>  Family: gaussian  ( identity )
#> Formula:          jaccard_scaled ~ 0 + year_nb/scaled_dist_up_km + (0 + year_nb |  
#>     main_bas/siteid) + (0 + year_nb | span) + (0 + year_nb:scaled_dist_up_km |      main_bas)
>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432
#> Data: data
#>  Offset: offset
#> 
#>      AIC      BIC   logLik deviance df.resid 
#>   1884.6   1946.9   -935.3   1870.6    54619 
#> 
#> Random effects:
#> 
#> Conditional model:
#>  Groups          Name                      Variance  Std.Dev.
#>  siteid.main_bas year_nb                   1.581e-04 0.012575
#>  main_bas        year_nb                   4.827e-05 0.006948
#>  span            year_nb                   1.055e-04 0.010271
#>  main_bas.1      year_nb:scaled_dist_up_km 8.619e-05 0.009284
#>  Residual                                  5.210e-02 0.228258
#> Number of obs: 54626, groups:  siteid:main_bas, 4916; main_bas, 291; span, 46
#> 
#> Dispersion estimate for gaussian family (sigma^2): 0.0521 
#> 
#> Conditional model:
<<<<<<< HEAD
#>                           Estimate Std. Error z value Pr(>|z|)    
#> year_nb                   0.026326   0.001795  14.669  < 2e-16 ***
#> year_nb:scaled_dist_up_km 0.007737   0.001651   4.687 2.77e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ggpredict(gaussian_jaccard_tmb$mod[[1]], "year_nb")
#> Warning: Following potential variables could not be found in the data: offset
#> Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 0 + year_nb | siteid:main_bas
#> Disable this warning with 'allow.new.levels=TRUE'
#> Warning: Following potential variables could not be found in the data: offset
#> # Predicted values of jaccard_dis_scaled
#> 
#> year_nb | Predicted |       95% CI
#> ----------------------------------
#>       0 |      0.00 | [0.00, 0.00]
#>      10 |      0.26 | [0.23, 0.30]
#>      20 |      0.53 | [0.46, 0.60]
#>      30 |      0.79 | [0.68, 0.90]
#>      40 |      1.05 | [0.91, 1.19]
#>      50 |      1.32 | [1.14, 1.49]
#>      60 |      1.58 | [1.37, 1.79]
#>      70 |      1.84 | [1.60, 2.09]
#> 
#> Adjusted for:
#> * jaccard_dis_scaled =  0.35
#> *  scaled_dist_up_km = -0.00
#> *             siteid = NA (population-level)
#> *           main_bas = NA (population-level)
#> *               span = NA (population-level)

======= #> Estimate Std. Error z value Pr(>|z|) #> year_nb -0.026326 0.001795 -14.669 < 2e-16 *** #> year_nb:scaled_dist_up_km -0.007737 0.001651 -4.687 2.77e-06 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ggpredict(gaussian_jaccard_tmb$mod[[1]], "year_nb") #> # Predicted values of jaccard_scaled #> # x = year_nb #> #> x | Predicted | 95% CI #> ------------------------------- #> 0 | 1.00 | [ 1.00, 1.00] #> 10 | 0.74 | [ 0.70, 0.77] #> 20 | 0.47 | [ 0.40, 0.54] #> 30 | 0.21 | [ 0.10, 0.32] #> 40 | -0.05 | [-0.19, 0.09] #> 50 | -0.32 | [-0.49, -0.14] #> 60 | -0.58 | [-0.79, -0.37] #> 70 | -0.84 | [-1.09, -0.60] #> #> Adjusted for: #> * scaled_dist_up_km = -0.00 #> * siteid = NA (population-level) #> * main_bas = NA (population-level) #> * span = NA (population-level) #> * offset = 1.00 >>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

rigal_trends_df %>%
  ggplot(aes(x = jaccard , y = hillebrand)) +
  geom_point() +
#  geom_smooth(method = "gam") +
  labs(x = "Jaccard trends (similarity, binary)", y = "Hillebrand trends
    (similarity, relative abundance)")

<<<<<<< HEAD

=======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

rigal_trends_df_loc <- filtered_dataset$location %>%
  left_join(rigal_trends_df, by = "siteid") %>%
  st_as_sf(coords = c("longitude", "latitude"),
  crs = 4326)

5.3 Relationships between temporal trends

knitr::include_graphics(here("doc", "fig", "p_cor_slope_tot.png"))

ti <- expand.grid(
  resp1 = unique(slope$response),
  resp2 = unique(slope$response)
  ) %>%
  filter(resp2 != resp1) %>%
  filter(
    resp1 %in% c("chao_richness", "species_nb", "log_species_nb",
      "total_abundance", "log_total_abundance")) %>%
  mutate_all(as.character) %>%
  arrange(resp1)

test <- map2(ti$resp1, ti$resp2,
  function(x, y) {
    bi <- slope %>%
      filter(response %in% c(x, y)) %>%
      select(siteid, response, linear_slope) %>%
      pivot_wider(names_from = "response", values_from = "linear_slope")

    return(bi)
  }
)

p_trends_trends <- map(test, function(x) {
  l <- colnames(x)
  x %>%
    ggplot(aes(x = !!sym(l[2]), y = !!sym(l[3]))) +
    geom_point() +
    geom_smooth(method = "loess") +
    labs(
      x = get_var_replacement()[l[2]],
      y = get_var_replacement()[l[3]]
    )
})
names(p_trends_trends) <- map_chr(test, ~colnames(.x)[2])

plot_grid(
  plotlist = p_trends_trends[names(p_trends_trends) %in% "log_species_nb"],
  ncol = 3
)
<<<<<<< HEAD
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

======= #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x' #> `geom_smooth()` using formula 'y ~ x'

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

5.4 Spatial trends

5.5 Rarity

Cf other document

5.6 Spatial scale

Selection on basin data

<<<<<<< HEAD

=======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

5.7 Environment

5.7.1 PCA

library(ade4)
library(factoextra)
#> Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
tar_load(rigal_slp_df)
slope_df <- rigal_slp_df

res.pca <- dudi.pca(select(slope_df, -siteid),
  scannf = FALSE,   # Hide scree plot
  nf = 5            # Number of components kept in the results
)

<<<<<<< HEAD

fviz_eig(res.pca)

=======

fviz_eig(res.pca)
#> Registered S3 methods overwritten by 'car':
#>   method                          from
#>   influence.merMod                lme4
#>   cooks.distance.influence.merMod lme4
#>   dfbeta.influence.merMod         lme4
#>   dfbetas.influence.merMod        lme4

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

p_pca <- map(list(c(1,2), c(2, 3), c(1, 3)),
  ~fviz_pca_var(res.pca,
    axes = .x,
    col.var = "contrib", # Color by contributions to the PC
    gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
    repel = TRUE     # Avoid text overlapping
    ) +
  theme(legend.position = "none") + labs(title = "")
)

plot_grid(plotlist = p_pca, ncol = 1)
<<<<<<< HEAD
#> Warning: ggrepel: 7 unlabeled data points (too many overlaps). Consider
#> increasing max.overlaps

======= #> Warning: ggrepel: 7 unlabeled data points (too many overlaps). Consider increasing max.overlaps

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

5.7.1.1 Environment

riv_data_pca <- riv %>%
  select(-siteid) %>%
  mutate(across(c(dist_up_km, dis_m3_pyr, urb_pc_cse, ele_mt_cav, pac_pc_cse),
      ~log(.x + 2))
  ) %>%
  rename(get_river_atlas_significant_var()) %>%
  na.omit()
pca_riv <- dudi.pca(df = scale(riv_data_pca), scannf = FALSE, nf = 3)

We see two main axis (as I found in France):
- First axis related to stream size (distance from source, River area)
- Second axis related to Temperature

fviz_eig(pca_riv)

<<<<<<< HEAD

=======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432


p_pca <- map(list(c(1,2), c(2, 3), c(1, 3)),
  ~fviz_pca_var(pca_riv,
    axes = .x,
    col.var = "contrib", # Color by contributions to the PC
    gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
    repel = TRUE     # Avoid text overlapping
    ) +
  theme(legend.position = "none") + labs(title = "")
)

p_pca[[1]]

<<<<<<< HEAD

p_pca[[2]]

p_pca[[3]]

=======

p_pca[[2]]

p_pca[[3]]

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432 #### Environment - Trends

tar_load(slp_env)
slp_env_for_pca <- slp_env %>%
  select(-all_of(c("siteid", "main_bas", "ecoregion", "geometry"))) %>%
  rename(get_river_atlas_significant_var())
slp_env_for_pca$geometry <- NULL
  
pca_slp_env <- dudi.pca(
  df = scale(na.omit(slp_env_for_pca)),
  scannf = FALSE, nf = 3)

p_slp_env_pca <- map(list(c(1,2), c(2, 3), c(1, 3)),

  ~fviz_pca_var(pca_slp_env,
    axes = .x,
    col.var = "contrib", # Color by contributions to the PC
    gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
    repel = TRUE     # Avoid text overlapping
    ) +
  theme(legend.position = "none") + labs(title = "")
)

p_slp_env_pca[[1]]
<<<<<<< HEAD
#> Warning: ggrepel: 10 unlabeled data points (too many overlaps). Consider
#> increasing max.overlaps

p_slp_env_pca[[2]]
#> Warning: ggrepel: 10 unlabeled data points (too many overlaps). Consider
#> increasing max.overlaps

p_slp_env_pca[[3]]
#> Warning: ggrepel: 11 unlabeled data points (too many overlaps). Consider
#> increasing max.overlaps

5.8 Model

======= #> Warning: ggrepel: 10 unlabeled data points (too many overlaps). Consider increasing max.overlaps

p_slp_env_pca[[2]]
#> Warning: ggrepel: 10 unlabeled data points (too many overlaps). Consider increasing max.overlaps

p_slp_env_pca[[3]]
#> Warning: ggrepel: 11 unlabeled data points (too many overlaps). Consider increasing max.overlaps

5.8 Model

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

5.8.1 Jaccard

tar_load(c(gaussian_jaccard_tmb, beta_jaccard_tmb))

tidy_gaus_jy  <- gaussian_jaccard_tmb %>%
  filter(year_var == "year_nb") %>%
  .$mod
names(tidy_gaus_jy) <- gaussian_jaccard_tmb[
  gaussian_jaccard_tmb$year_var == "year_nb",
  ]$var_jaccard
dw <- dwplot(tidy_gaus_jy,
  vline = geom_vline(xintercept = 0, colour = "grey50", linetype = 2)) +
  xlab("Coefficient")
dw

<<<<<<< HEAD

=======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

diagnose(gaussian_jaccard_tmb$mod[[1]])
#> model looks OK!
diagnose(gaussian_jaccard_tmb$mod[[3]])
#> model looks OK!
diagnose(gaussian_jaccard_tmb$mod[[5]])
<<<<<<< HEAD
#> model looks OK!

pred <- ggemmeans(gaussian_jaccard_tmb$mod[[1]],
  terms = c("year_nb", "scaled_dist_up_km [quart2]"),
  type = "re")
#> Warning: Following potential variables could not be found in the data: offset
#> NOTE: A nesting structure was detected in the fitted model:
#>     scaled_dist_up_km %in% year_nb
#> Warning: Following potential variables could not be found in the data: offset
plot(pred)

======= #> Unusually large Z-statistics (|x|>5): #> #> year_nb:scaled_dist_up_km #> -5.182309 #> #> Large Z-statistics (estimate/std err) suggest a failure of the Wald approximation - often also associated with parameters that are at #> or near the edge of their range (e.g. random-effects standard deviations approaching 0). While the Wald p-values and standard errors #> listed in summary() are unreliable, profile confidence intervals (see ?confint.glmmTMB) and likelihood ratio test p-values derived by #> comparing models (e.g. ?drop1) may still be OK. (Note that the LRT is conservative when the null value is on the boundary, e.g. a #> variance or zero-inflation value of 0 (Self and Liang 1987; Stram and Lee 1994; Goldman and Whelan 2000); in simple cases the p-value #> is approximately twice as large as it should be.)

pred <- ggemmeans(gaussian_jaccard_tmb$mod[[1]],
  terms = c("year_nb", "scaled_dist_up_km [quart2]"),
  type = "re")
#> NOTE: A nesting structure was detected in the fitted model:
#>     scaled_dist_up_km %in% year_nb
plot(pred)

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432


pred <- ggemmeans(gaussian_jaccard_tmb$mod[[3]],
  terms = c("year_nb", "scaled_dist_up_km [quart2]"),
  type = "re.zi")
<<<<<<< HEAD
#> Warning: Following potential variables could not be found in the data: offset
#> NOTE: A nesting structure was detected in the fitted model:
#>     scaled_dist_up_km %in% year_nb
#> Warning: Following potential variables could not be found in the data: offset
plot(pred)

======= #> NOTE: A nesting structure was detected in the fitted model: #> scaled_dist_up_km %in% year_nb plot(pred)

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432


pred <- ggemmeans(gaussian_jaccard_tmb$mod[[5]],
  terms = c("year_nb", "scaled_dist_up_km [quart2]"),
  type = "re.zi")
<<<<<<< HEAD
#> Warning: Following potential variables could not be found in the data: offset
#> NOTE: A nesting structure was detected in the fitted model:
#>     scaled_dist_up_km %in% year_nb
#> Warning: Following potential variables could not be found in the data: offset
plot(pred)

======= #> NOTE: A nesting structure was detected in the fitted model: #> scaled_dist_up_km %in% year_nb #> Warning: Can't compute random effect variances. Some variance components equal zero. Your model may suffer from singularity (see '?lme4::isSingular' and #> '?performance::check_singularity'). #> Solution: Respecify random structure! You may also decrease the 'tolerance' level to enforce the calculation of random effect variances. plot(pred)

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

Map:

hist(coef(gaussian_jaccard_tmb$mod[[1]])$cond[["siteid:main_bas"]]$year_nb)

<<<<<<< HEAD

hist(coef(gaussian_jaccard_tmb$mod[[1]])$cond[["main_bas"]]$year_nb)

=======

hist(coef(gaussian_jaccard_tmb$mod[[1]])$cond[["main_bas"]]$year_nb)

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

tar_load(filtered_dataset)
loc <- filtered_dataset$location %>%
  st_as_sf(coords = c("longitude", "latitude"),
  crs = 4326)
world <- ne_countries(scale = "medium", returnclass = "sf")

pred_slope <- gaussian_jaccard_tmb %>%
  mutate(
    site_slope = map(mod, function(m) {
      coef(m)$cond[["siteid:main_bas"]] %>%
        rownames_to_column(var = "siteid_main_bas") %>%
        separate(
          col = "siteid_main_bas",
          into = c("siteid", "main_bas"),
          sep = ":") %>%
        select(-main_bas)
    }),
    basin_slope = map(mod, function(m) {
      coef(m)$cond[["main_bas"]] %>%
        rownames_to_column(var = "main_bas") 
    }),
  ) %>%
select(-mod)
pred_slope$basin_slope[[1]] %>%
  as_tibble()
<<<<<<< HEAD
#> # A tibble: 291 x 3
#>    main_bas   year_nb `year_nb:scaled_dist_up_km`
#>    <chr>        <dbl>                       <dbl>
#>  1 1080023890  0.0288                     0.00504
#>  2 1080040200  0.0260                     0.00331
#>  3 2080008490  0.0441                     0.00582
#>  4 2080016240  0.0198                     0.0113 
#>  5 2080016250  0.0280                     0.00676
#>  6 2080016280  0.0202                     0.00793
#>  7 2080016290  0.0170                     0.0120 
#>  8 2080016350  0.0306                     0.00758
#>  9 2080016400  0.0225                     0.00945
#> 10 2080016460  0.0292                     0.00693
#> # ... with 281 more rows
=======
#> # A tibble: 291 × 3
#>    main_bas   year_nb `year_nb:scaled_dist_up_km`
#>    <chr>        <dbl>                       <dbl>
#>  1 1080023890 -0.0288                    -0.00504
#>  2 1080040200 -0.0260                    -0.00331
#>  3 2080008490 -0.0441                    -0.00582
#>  4 2080016240 -0.0198                    -0.0113 
#>  5 2080016250 -0.0280                    -0.00676
#>  6 2080016280 -0.0202                    -0.00793
#>  7 2080016290 -0.0170                    -0.0120 
#>  8 2080016350 -0.0306                    -0.00758
#>  9 2080016400 -0.0225                    -0.00945
#> 10 2080016460 -0.0292                    -0.00693
#> # … with 281 more rows
>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

pred_slope <- pred_slope %>%
  mutate(
    sf_site = map(site_slope,
      ~left_join(.x, select(loc, siteid, geometry), by = "siteid") %>%
        st_as_sf()
      ),
    sf_basin = map(basin_slope,
      ~left_join(.x, loc %>%
        mutate(main_bas = as.character(main_bas)) %>%
        select(main_bas, geometry), by = "main_bas") %>%
        st_as_sf()
      ),
    sf_site_map = map2(sf_site, year_var,
      ~ggplot(data = world) +
        geom_sf() +
        geom_sf(data = .x, aes(color = year_nb), shape = 16) +
        scale_color_viridis() +
        theme(legend.position = "bottom")),
    sf_basin_map = map2(sf_basin, year_var,
      ~ggplot(data = world) +
        geom_sf() +
        geom_sf(data = .x, aes(color = year_nb), shape = 16) +
        scale_color_viridis() +
        theme(legend.position = "bottom"))
  )

<<<<<<< HEAD

5.8.1.1 Site

pred_slope$sf_site_map[[1]] +
  labs(title = "Jaccard")


pred_slope$sf_site_map[[3]] +
  labs(title = "Turnover")


pred_slope$sf_site_map[[5]] +
  labs(title = "Nestedness")

5.8.1.2 Basin

pred_slope$sf_basin_map[[1]] +
  labs(title = "Jaccard")


pred_slope$sf_basin_map[[3]] +
  labs(title = "Turnover")


pred_slope$sf_basin_map[[5]] +
  labs(title = "Nestedness")

=======

5.8.1.1 Site

pred_slope$sf_site_map[[1]] +
  labs(title = "Jaccard")


pred_slope$sf_site_map[[3]] +
  labs(title = "Turnover")


pred_slope$sf_site_map[[5]] +
  labs(title = "Nestedness")

5.8.1.2 Basin

pred_slope$sf_basin_map[[1]] +
  labs(title = "Jaccard")


pred_slope$sf_basin_map[[3]] +
  labs(title = "Turnover")


pred_slope$sf_basin_map[[5]] +
  labs(title = "Nestedness")

5.9 Species richness

tar_load(gaussian_rich_tmb)

lry <- gaussian_rich_tmb$mod[[5]]
cry <- gaussian_rich_tmb$mod[[1]]

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

knitr::opts_chunk$set(eval = FALSE)

paste0(var_temporal_trends,
  " ~ ",
  "dist_up_km + tmp_dc_cyr +
  (1 + dist_up_km + tmp_dc_cyr | main_bas)"
) %>%
as.formula

library(spaMM)
tar_load(c(trend_env))
names(trend_env) <- var_temporal_trends

ti_trends <- map(var_temporal_trends,
  ~fitme(
  as.formula(paste0("scale(", .x, ") ~ dist_up_km + tmp_c_cyr + (1 | main_bas)")),
  data = slp_env[,
    colnames(slp_env) %in%
    c(.x, "siteid", "dist_up_km", "tmp_c_cyr", "main_bas", "ecoregion")
  ]
  )
)
names(ti_trends) <- var_temporal_trends

ci <- map_dfr(ti_trends[names(ti_trends) %in% c("log_total_abundance", "log_species_nb", "jaccard",
      "appearance", "disappearance", "turnover",
      "nestedness", "hillebrand")], function(v) {
  x <- confint(v, c("dist_up_km", "tmp_c_cyr"), verbose = FALSE)
  tmp <- map_dfr(x, function(x) {
    tm <- x$interval
    names(tm) <- c("lower", "upper")
    return(tm)
  }, .id = "parm"
  )
  return(tmp)
}, .id = "response")
ci %>%
  kable(.,
    caption = "Confidence interval (bootstrap)")

from_data_to_predict <- function(
  model_list = NULL,
  dataset = NULL,
  re = NULL
) {

  pred <- na.omit(dataset)
  for (i in seq_along(names(model_list))) {
    tmp <- predict(model_list[[i]],
      re.form = as.formula(re),
      intervals = "predVar"
      )[,1]
    names(tmp) <- NULL
    pred[[names(model_list)[i]]] <- tmp
  }
  return(pred)
}

xx <- list(
    tmp_c_cyr = seq(min(na.omit(slp_env$tmp_c_cyr)), max(na.omit(slp_env$tmp_c_cyr)), length.out = 10),
    dist_up_km = seq(min(na.omit(slp_env$dist_up_km)), max(na.omit(slp_env$dist_up_km)), length.out = 10),
    main_bas = unique(na.omit(slp_env$main_bas))
)
new_data <- expand.grid(xx) %>%
    as_tibble()
pred <- predict(ti_trends[[1]], newdata = new_data, re.form = re, intervals="predVar")
int <- cbind(pred[,1], attr(pred, "intervals")) %>%
    as_tibble() %>%
    rename(y = V1)

xxx <- cbind(int, new_data) %>%
    as_tibble
xxx %>%
    ggplot(aes(y = y, x = dist_up_km, color = main_bas)) +
    geom_point() +
    theme(legend.position = "none")

re <- paste0("~ tmp_c_cyr + (1 | main_bas)")
pred_tmp <- from_data_to_predict(
  model_list = ti_trends,
  dataset = slp_env,
  re = re
)
p_tmp <- map(
  c(
      "log_total_abundance", "log_species_nb", "jaccard",
      "appearance", "disappearance", "turnover",
      "nestedness", "hillebrand"),
  function(x) {
    pred_tmp %>%
      ggplot(aes(
          y = !!sym(x),
          x = tmp_c_cyr,
          color = main_bas
          )) +
      geom_line() +
      theme(legend.position = "none")
  }
)
plot_grid(plotlist = p_tmp, ncol = 3)

re <- paste0("~ dist_up_km + (1 | main_bas)")
pred <- na.omit(slp_env)
for (i in seq_along(var_temporal_trends)) {
  tmp <- predict(trend_env[[i]],
    re.form = as.formula(re)
    )[, 1]
  names(tmp) <- NULL
  pred[[var_temporal_trends[i]]] <- tmp
}
p <- map(
  c("log_total_abundance", "log_species_nb", "jaccard",
      "appearance", "disappearance", "turnover",
      "nestedness", "hillebrand"),
  function(x) {

    pred %>%
      ggplot(aes(
          y = !!sym(x),
          x = dist_up_km,
          color = as.factor(main_bas),
          shape = as.factor(ecoregion)
          )) +
      geom_line() +
      theme(legend.position = "none")
  }
)
plot_grid(plotlist = p, ncol = 3)

re <- paste0("~ tmp_dc_cyr + (1 | main_bas)")
pred_tmp <- from_data_to_predict(
  model_list = trend_env,
  dataset = slp_env,
  re = re
)
p_tmp <- map(
  c("log_species_nb", "jaccard", "appearance", "disappearance", "turnover", "hillebrand"),
  function(x) {
    pred_tmp %>%
      ggplot(aes(
          y = !!sym(x),
          x = tmp_dc_cyr,
          color = as.factor(main_bas)
          )) +
      geom_line() +
      theme(legend.position = "none")
  }
)
plot_grid(plotlist = p_tmp, ncol = 3)

Ecoregion:

re <- paste0("~ dist_up_km + (1 + dist_up_km | ecoregion)")
pred_dist_ecoregion <- from_data_to_predict(
  model_list = trend_env,
  dataset = slp_env,
  re = re
)

p_dist_ecoregion <- map(
  c("log_species_nb", "jaccard", "appearance", "disappearance", "turnover", "hillebrand"),
  function(x) {

    pred_dist_ecoregion %>%
      ggplot(aes(
          y = !!sym(x),
          x = dist_up_km,
          color = as.factor(ecoregion)
          )) +
      geom_line() +
      theme(legend.position = "none")
  }
)
plot_grid(plotlist = p_dist_ecoregion, ncol = 3)

re <- paste0("~ tmp_dc_cyr + (1 + tmp_dc_cyr | ecoregion)")
pred_tmp_ecoregion <- from_data_to_predict(
  model_list = trend_env,
  dataset = slp_env,
  re = re
)

p_tmp_ecoregion <- map(
  c("log_species_nb", "jaccard", "appearance", "disappearance", "turnover", "hillebrand"),
  function(x) {
    pred_tmp_ecoregion %>%
      ggplot(aes(
          y = !!sym(x),
          x = tmp_dc_cyr,
          color = as.factor(main_bas)
          )) +
      geom_line() +
      theme(legend.position = "none")
  }
)
plot_grid(plotlist = p_tmp_ecoregion, ncol = 3)

<<<<<<< HEAD

5.9 Analysis

5.10 Reproducibility

=======

5.10 Analysis

5.11 Reproducibility

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

Reproducibility receipt

## datetime
Sys.time()

## repository
if(requireNamespace('git2r', quietly = TRUE)) {
  git2r::repository()
} else {
  c(
    system2("git", args = c("log", "--name-status", "-1"), stdout = TRUE),
    system2("git", args = c("remote", "-v"), stdout = TRUE)
  )
}

## session info
sessionInfo()

References

<<<<<<< HEAD

Blowes, Shane A., Sarah R. Supp, Laura H. Antão, Amanda Bates, Helge Bruelheide, Jonathan M. Chase, Faye Moyes, et al. 2019. “The Geography of Biodiversity Change in Marine and Terrestrial Assemblages.” Science 366 (6463): 339–45. https://doi.org/10.1126/science.aaw1620.

Dornelas, Maria, Nicholas J. Gotelli, Brian McGill, Hideyasu Shimadzu, Faye Moyes, Caya Sievers, and Anne E. Magurran. 2014. “Assemblage Time Series Reveal Biodiversity Change but Not Systematic Loss.” Science 344 (6181): 296–99. https://doi.org/10.1126/science.1248484. =======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

Hillebrand, Helmut, Bernd Blasius, Elizabeth T. Borer, Jonathan M. Chase, John A. Downing, Britas Klemens Eriksson, Christopher T. Filstrup, et al. 2018. “Biodiversity Change Is Uncoupled from Species Richness Trends: Consequences for Conservation and Monitoring.” Journal of Applied Ecology 55 (1): 169–84. https://doi.org/https://doi.org/10.1111/1365-2664.12959.

Klink, Roel van, Diana E. Bowler, Konstantin B. Gongalsky, Ann B. Swengel, Alessandro Gentile, and Jonathan M. Chase. 2020. “Meta-Analysis Reveals Declines in Terrestrial but Increases in Freshwater Insect Abundances.” Science 368 (6489): 417–20. https://doi.org/10.1126/science.aax9931.

<<<<<<< HEAD

Vellend, Mark, Lander Baeten, Isla H. Myers-Smith, Sarah C. Elmendorf, Robin Beauséjour, Carissa D. Brown, Pieter De Frenne, Kris Verheyen, and Sonja Wipf. 2013. “Global Meta-Analysis Reveals No Net Change in Local-Scale Plant Biodiversity over Time.” Proceedings of the National Academy of Sciences 110 (48): 19456–59. =======

>>>>>>> bd3c5f9d6923adad81aa20a3de91333a12b80432

Meeting report

Alain Danet

11 February, 2022

10 février, 2022

1 Introduction

1.1 Why do we want to assess community changes ?

1.2 Documenting community changes

1.3 Research gap

1.4 Research questions

1.4 Research questions

1.4.1 Drivers of temporal changes

1.4.2 Within and basin scale

1.5 Fish freshwater in streams

2 Methods

2.1 Data selection

2.2 Biodiversity

2.2.1 Quantity

2.2.2 Betadiversity

3 Environment

3 Environment

3.1

4 Statistical analysis

4.1 Species composition change

4.2 Model

4.2.1 Turnover

4.2 Model

4.2.1 Turnover

5 Results

5.1 Direction of trends

5.2 Average

5.3 Relationships between temporal trends

5.4 Spatial trends

5.5 Rarity

5.6 Spatial scale

5.7 Environment

5.7.1 PCA

5.7.1.1 Environment

5.8 Model

5.8 Model

5.8.1 Jaccard

5.8.1.1 Site

5.8.1.2 Basin

5.8.1.1 Site

5.8.1.2 Basin

5.9 Species richness

5.9 Analysis

5.10 Reproducibility

5.10 Analysis

5.11 Reproducibility

References